%%  Skrypt dla identyfikacji systemu
% Parametry:
%   u       -   wejscie do systemu
%   y       -   wyjscie z systemu
%   Ts      -   czas probkowania
%   na      -   ilosc wezlow
%   nb      -   ilosc zer i 1
%   nk      -   opoznienia
%   nb i nk -   co do wymiarow takie samo
% Wyjscie:
%   sys     -   system liniowy reprezentujacy relacje pomiedzy u i y
%
%   sys:
%       x' = Ax + Bu
%       y = Cx + Du
%

function sys = identyfikuj(u,y,Ts,na,nb,nk)

%% przygotowanie danych
if length(u(:,1)) < 50
    u = u-ones(length(u(:,1)),1)*mean(u(1:end,:));
    y = y-ones(length(y(:,1)),1)*mean(y(1:end,:));
else
    u = u-ones(length(u(:,1)),1)*mean(u(1:50,:));
    y = y-ones(length(y(:,1)),1)*mean(y(1:50,:));
end

%% tworzenie obiektow danych dla identyfikacji

ze = iddata(y,u,Ts);
ny = length(y(1,:));
nu = length(u(1,:));

%% wyznaczenie modelu arx

orders = [na nb nk];
armaxmodel = arx(ze,orders);

%% wyznaczenie systemu

sys = ss(armaxmodel);
sys = sys(1:ny,1:nu);